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Abstract: 

The absence of national local geoid model in some countries has led to the determination of local geoid model in various parts 
of those countries. Local geoid models are determined using the geometric and gravimetric methods amongst others. Using 
the geometric method requires fitting an interpolation surface to points of known geoidal undulations which requires the 
determination of the geometric geoid model parameters and its accuracy using least squares technique. Because of the 
rigorous as well as the matrix nature of the technique, researcher have been experiencing difficulty in its application for the 
determination of geometric geoid models’ parameters and their accuracy. This paper presents a detailed procedure for the 
determination of geometric geoid models' parameters as well as their accuracy using least squares technique. The steps to be 
considered when applying the technique are enumerated in sequential order. The enumerated steps were also demonstrated 
with a numerical example. 

Keywords: Geometric geoid, model parameters, accuracy, least squares 


1. Introduction 

Least squares are a statistical method used to determine a line of best fit by minimizing the sum of squares created by 
a mathematical function. It is a popular method for determining regression equations. Instead of trying to solve an equation 
exactly, least squares method is used to determine a close approximation which is known as the estimate. Modelling methods 
that are often used when fitting a function to a curve include the straight-line method, polynomial method, logarithmic method 
and Gaussian method. The Least-Squares criterion is an imposed condition for obtaining a unique solution for an incompatible 
system of linear equations. The term adjustment, in a statistical sense, is a method of deriving estimates for random variables 
from their observed values. The application of the least-squares criterion in the adjustment problem is called the Least- 
Squares Adjustment method (Mohammad-Karim, 1981). The method of least squares is a rigorous technique that can be 
applied to the adjustment of horizontal geodetic network to yield the most likely values of the survey measurements. In 
geodesy, it is desirable or necessary to fit a plane or curve surface to a set of points with known coordinates or heights. In 
solving this type of problem, it is first necessary to decide on the appropriate functional model for the data as stated by Ghilani 
(2010). The decision as to whether to use a plane or curve surface depends on the size of the application area. To determine 
the best fit surface, two or more surfaces have to be applied and the one with smaller residuals after least squares solution 
with the surfaces selected. 

Geometric geoid models are surfaces that are fitted to the geoidal undulations of an area to enable geoid heights of 
new points within the area to be interpolated. These surfaces are plane as well as curve surfaces depends on the degree. The 
curve surfaces are ether quadratic or polynomial in nature. The plane surfaces are usually applied in small areas while the 
curve surfaces are applied in relatively large areas. The larger the area the higher the order as well as the degree of the 
polynomial model/surface. To apply any of these models in a particular area, the model parameters as well as its accuracy 
have to be determined using least squares technique. Obtaining the accuracy of the model enables the reliability of the model 
to be determined. 

Various researchers have been experiencing difficulty in the application of least squares adjustment technique for 
determination of geometric geoid models’ parameters and their accuracy. The difficulty in its application resulted from its 
matrix nature. The computation of these parameters cannot be handled by Least squares adjustment software as the model 
terms are not obtained directly from measurement, that is, they are not bearings, azimuth, angles, distances, change in 
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northing and easting. They are normally reduced from the variables. Also, the dimension of the coefficient matrix depends on 
the number of points heights and the order as well as the degree of the model. The general matrix notation of least squares 
adjustment is simply the sum of the estimate and the matrix of observations equals to the residual matrix. In geoid model 
parameters determination, the residual is not truly needed. Making the matrix notation of least squares adjustment for 
geometric geoid parameters determination to be simply observation equal the estimate. Previous studies in which least 
squares adjustment technique was applied for the determination of local geometric geoid models never presented the 
breakdown of how the technique was applied in the studies. To determine the reliability of the model, the geoid heights of 
points from the model are compared with their known geoid heights. 

This paper presents step by step application of observation equation method of least squares adjustment technique 
for determination of local geometric geoid models’ parameters and their accuracy. 

1.1. Observation Equation Method of Least Squares Adjustment 

Equations that relate observed quantities to both observational residuals and independent unknown parameters are 
called observation equations. One equation is written for each observation and for a unique set of unknowns. For a unique 
solution of unknowns, the number of equations must equal the number of unknowns. Usually, there are more observations 
(and hence equations) than unknowns, and this permits determination of the most probable values for the unknowns based on 
the principle of least squares (Ghilani and Wolf, 2006). 

Ayeni (2001) and Okwuashi and Asuquo (2014) explained that, in the observations equation method, the adjusted 
observations are expressed as a function of the adjusted parameter. The functional relationship between adjusted 
observations and the adjusted parameters as given in Ono et al (2014) is: 

L a = F(X a ) ( 1 ) 

Where, L a = adjusted observations and X a = adjusted parameters. Equation (1) is linear function and the general 
observation equation model was obtained. 

To make the matrix expression for performing least squares adjustment, analogy will be made with the systematic 
procedures. The system of observation equations is presented by matrix notation as (Mishima and Endo, 2002 and Ono et al, 
2018): 


V = AX - L 


( 2 ) 


But the residual matrix is not necessary when applying least squares adjustment technique for the determination of local 
geometric geoid model parameters. Thus, the general matrix notation becomes 


L = AX 

where, A = Design Matrix, X = Vector of Unknowns, L = Observation Matrix. That is, 


( 3 ) 
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The determination of the unknown parameters, X, requires the normal matrix, N and the matrix of numeric terms, t to be 
deduced. It is to be noted here that the observations are not weighted. 

According to Ghilani (2010), a system of unweighted linear observation equations can be expressed in matrix notation as: 

A t AX=A t L (4) 

To make X the subject of the formula, both sides of equation (4) will be divided by A r A. Thus, 

X = (A 7 A)~ l A T L (5) 

If A t A = N, normal matrix and A r L = t, matrix of numeric terms, then equation (5) becomes 
X = N^t (6) 

1.2. Accuracy/Reliability of Geometric Geoid Model 

The accuracy of determined local geometric geoid model is obtained using the Root Mean Squares Error, RMSE index. 
To evaluate the determined local geometric geoid model accuracy, the local geoid model is used to determine the geoidal 
heights of points whose geoid heights are known. The geometric geoid model geoidal undulations are compared with the 
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known geoidal undulations of the points to obtain the residuals. The Root Mean Squares Error, RMSE index for the 
computation of geometric geoid model accuracy as given by Kao et al (2017} is 



( 7 } 


Where, 


V — Nknown Nmodel (Residual) 

Nknown = Known geoid height of point 
Nmodel = Model geoid height of point 
n = Number of points 

1.3. Bicubic Geoid Model 

The bicubic surface is one of the polynomial geometric geoid models and it is given as 

N = a o + ci\X + a 2 y + a 3 x 2 + a 4 y 2 + a 5 xy + a 6 x 2 y + a 1 xy 2 + a s x 3 + a 9 y 3 , (8) 


Where, 


y = ABS(Y-Y a ) [g] 

x = ABS(X - X a ) 

Y = Northing coordinates of observed station 
X = Easting coordinates of observed station 

Y o = Northing coordinates of the origin (average of the northing coordinates) 

X 0 = Easting coordinates of the origin (average of the easting coordinates) 

2. Steps to Be Considered When Computing Geometric Geoid Models Parameters and Accuracy Using Observation 
Equation Method of Least Squares Adjustment Technique 

The steps to be considered when computing geometric geoid models’ parameters and their accuracy using the least 
squares method are as follow: 

• Deduce the coefficient matrix, A, observation matrix, L and matrix of unknown parameters, X from the given or model 
of interest. It is to be noted here that the geometric geoid model is already established. The coefficients of the model 
terms on one side of the equation form the matrix of unknown parameters while the terms form the coefficient matrix. 
The model terms which form the coefficient matrix are deduced from the given variables. The observation matrix is 
deduced from the term on the other side of the equation as well as the model. The number of points heights deduced 
with respect to the centroid of the study area determines the number of observation equations. In this, the weights of 
the observed heights are not truly necessary. 

• Having deduced the above stated matrices, the model parameters are computed using equations (6). 

• Since the parameters are computed, the next step is to substitute the computed parameters accordingly in the model. 

• Having substituted the computed parameters in the model, it (the model) can now be used to develop a 

program such that the geoid heights of new points in the study area can be determined if the variables of the new 
points are given. 

• The model accuracy has to be computed using equation (7) to determine its reliability. To determine the model 
accuracy, the geoidal heights of points obtained from the model are compared with their corresponding known geoid 
heights to obtained the residuals. 

2.1. Numerical Application in Bicubic Geoid Model Parameters and Accuracy Determination 

The ellipsoidal heights of eight benchmarks whose rectangular coordinates and orthometric heights are known were 
determined using DGPS. The rectangular coordinates, orthometric heights, ellipsoidal heights and the geoid heights as 
computed from the orthometric and ellipsoidal heights are given in table 1. To interpolate geoid heights of new points, the 
bicubic geoid model, equation (8) is to be fitted to the computed geoid heights of the benchmarks. Compute the model 
parameters and accuracy using least squares technique. 
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Station 

Coordinates 

Orthometric 
Height, H 

Ellipsoidal 
Height, h 

Geoidal Height, 

N (m) 

Northing 

Easting 

ESOOl 

249308.287 

354033.425 

175.189 

209.237 

34.048 

ESO02 

244533.051 

278026.486 

291.686 

326.581 

34.895 

ESO03 

249241.822 

362785.077 

222.300 

256.677 

34.377 

ESO04 

259174.974 

355889.303 

425.449 

460.033 

34.584 

ESO05 

247210.935 

359597.719 

325.386 

359.665 

34.279 

ESO06 

260606.174 

332700.238 

120.829 

155.115 

34.286 

ESO07 

252751.094 

344865.087 

143.546 

177.828 

34.282 

ESO08 

276864.558 

374129.027 

315.314 

349.689 

34.375 

ESO09 

243587.154 

340245.247 

257.359 

291.966 

34.607 

ESO 10 

269356.441 

361478.369 

199.075 

233.531 

34.456 

ESO 11 

256457.248 

357864.254 

351.273 

385.448 

34.175 


Table 1: Coordinates and Heights of Benchmarks 


Solution 

Using equation (9), Y o and X o are computed as the mean as well as the centroid of the given positions as: 
Y a =255371.976 
X 0 =347419.476 

Deduction of normal equations using equation (8] 


2 2 2 2 3 3 

N x - a g + a x x x + a 2 y x + a 3 x x + a 4 y{ + a 5 xy x + a 6 x y x + a 1 xy 1 + a 8 x x + a 9 ^ 


N 2 = ci 0 + a x x 2 + a 2 y 2 + a 3 xl + a 4 y 2 + a 5 xy 2 + a 6 x 2 y 2 + a 1 xy\ + a s x \ + a 9 y\ 

2 2 2 2 3 3 

N 3 - a o + a x x 3 + a 2 y 3 + ci 3 x 3 + a 4 y 3 + a 5 xy 3 + a 6 x~y 3 + a 1 xy 3 + a & x 3 + a 9 y 3 
N 4 =a 0 + a x x 4 +a 2 y 4 +a 3 x 2 +a 4 y 4 +a 5 xy 4 +a 6 x 2 y 4 + a 7 xy 2 4 +a 8 x 4 +a 9 y 4 
N 5 =a 0 + a { x 5 +a 2 y 5 + a 3 x] + a 4 yl + a 5 xy 5 +a 6 x 2 y 5 + a 7 xy 2 +a 8 x 3 + a 9 y] 

N 6 =a 0 + a x x 6 + a 2 y 6 + a 3 x 2 + a 4 y\ + a 5 xy 6 + a 6 x 2 y 6 + a 7 xy 2 + a 8 x 6 3 + a 9 y\ 

N 7 = a g + a x x 7 + a 2 y 7 + a 3 x 2 + a 4 y 2 + a 5 xy 7 + a 6 x 1 y 1 + a 7 xy 7 + a s x 2 + a 9 y 2 
A/g = a D + a x x s +a 2 y 8 + a 3 x % + a 4 y s + a 5 xy s + a 6 x~y % + a 7 xy 8 + a 8 x 8 + a 9 y % 

N 9 - a n + a x x 9 + a 2 y 9 + a 3 x 9 + a 4 y 9 + a 5 xy 9 + a 6 x y 9 + a 1 xy g + a s x 9 + a 9 y g 
N w = a 0 + a x x 10 + a 2 y l0 + a 3 xf 0 + a 4 y 2 w + a 5 xy l0 + a b x 2 y m + a 7 xy 2 0 + a 8 xf 0 + a g y x0 
N n = a 0 + a x x n + a 2 y n + a 3 xf x + a 4 yl x d-aj-ryn +a b x 2 y n + a 7 xyl x + a 8 xj\ +a g y xl 


Also using equation (3), the coefficient matrix, A, matrix of unknown parameters, X and observation matrix, L are: 


A = 



y i 

A 

y i 

*>i 

x 2 y x 

xy 2 

A 

>i 3 

X, 

y 2 

X 2 

y 2 

xy 2 

x 2 y 2 

xv 2 

x\ 

yl 

x 3 

y 3 

x\ 

>3 

.ry 3 

x 2 y 3 

xy 2 

x\ 

y 1 

X 4 

y 4 

A 

> 4 

*y 4 

x 2 y 4 

xyl 

A 

y l 

x 5 

y s 

A 

>5 

xy 5 

•* 2 >s 

xy 2 

A 

yl 

*6 

y 6 

A 

>6 

xy 6 

* 2 .v 6 

xyl 

A 

yl 

*7 

y? 

A 

> 7 

xy 7 

x 2 y 7 

xy 2 

A 

yl 

*8 

y 8 

A 

> 8 2 

xy s 

x 2 y g 

xyl 

xl 

yl 

x 9 

>9 

Xg 

>9 

xy g 

x 2 y g 

xyl 

A 

yl 

*10 

>10 

Ao 

>10 

xy 10 

* 2 y io 

xy 2 w 

v 3 

A 10 

>10 

x n 

>11 

Ai 

>U 

-*>11 

-* 2 .Vii 

xyh 

An 

>n 


L = 


N, 

N 2 

n 3 

n 4 

n 5 

n 6 

N, 

N, 

N g 

N l0 


f 34 . 048 ^ 

34.895 

34.377 

34.584 

34.279 

34.286 

34.282 

34.375 

34.607 

34.456 

34.175 


f a ^ 


X = 


\ U 9j 
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Considering equation (3), the above matrices are rearranged as 


L = A 


aO 


'l 


yi 

2 

*1 

2 

yi 

xy i 

2 

X 

2 

XJl 

x i 

y 2 

n 2 


1 

*2 

y 2 

2 

x 2 

y 2 

*y 2 

* 2 y 2 

xy 2 

x\ 

y 2 

N 3 


1 

JC 3 

y 3 

Xj 

y 2 

xy 3 

x 2 .v 3 

xy 2 

X 3 

y 3 



1 

*4 

y 4 

*4 2 

yl 

xy 4 

* 2 y 4 

xy 4 


y 4 3 

n 5 


1 

*5 

?5 

* 5 2 

y 5 

xy 5 

x 2 y 5 

xy 2 

* 5 3 

yl 

n 6 

= 

1 

X 6 

y 6 

*1 

y 6 

*y 6 

x 2 y 6 

xy 6 

X 3 

yl 

n 7 


1 

x 7 

y 2 

2 

x 7 

2 

y 2 

xy 7 

2 

x y 7 

2 

xy 7 

x 7 

yl 

Ns 


1 

x & 

y 8 

2 

*8 

2 

^8 

xy* 

2 

x y 8 

2 

*y 8 

xl 

yl 

n 9 


1 

Xg 


2 

x 9 

yl 

xy 9 

2 

x y 9 

2 

xy 9 

x 9 

yl 

A^io 


1 

*10 

Vio 

X 2 

-'lO 

yi 

xy w 

x 2 y w 

xy 2 0 

X 3 

A 10 

yl 0 

N n j 



x n 

yu 

X 2 

-'ll 

yi 


x 2 y, i 

xyi 

An 

yi 


x X. 


f a f 


a , 


x 


a. 


a A 


a. 


a n 


\ a 9j 


Deduction of values of x and y for each equation using equation (9) and the computation of coefficient matrix, A elements as 
presented in table 2. 


S/N 

3o 

X 

y 

X 2 

y2 

xy 

x 2 y 

xy 2 

X 3 

y3 

1 

1 

6613. 

949 

6063. 

689 

43744326 

.185 

367683 

26.494 

401049 

33.205 

265251997 

452.369 

243183849 

614.972 

28932275 

8332.3 

222951703 

593.398 

2 

1 

69392 

.990 

10838 

.925 

48153870 

10.672 

117482 

299.097 

752145 

422.811 

521936195 

30177.900 

815244796 

3697.900 

33415410 

0926675.0 

127338185 

0100.930 

3 

1 

15365 

.601 

6130. 

154 

23610170 

5.266 

375787 

90.293 

941935 

05.455 

144733985 

5871.900 

577420711 

367.821 

36278446 

84395.0 

230363778 

461.463 

4 

1 

8469. 

827 

3802. 

998 

71737975 

.570 

144627 

92.405 

322107 

34.984 

272819364 

572.745 

122497354 

867.272 

60760826 

8493.0 

550019679 

61.395 

5 

1 

12178 

.243 

8161. 

041 

14830961 

1.424 

666025 

93.171 

993871 

45.613 

121036084 

6490.330 

811102588 

289.749 

18061505 

41087.3 

543546505 

687.124 

6 

1 

14719 

.238 

5234. 

198 

21665595 

6.596 

273968 

26.800 

770434 

01.522 

113402013 

5309.530 

403260404 

149.425 

31890105 

10466.4 

143400411 

060.929 

7 

1 

2554. 

389 

2620. 

882 

6524901. 

306 

686902 

3.411 

669475 

1.662 

171009975 

69.926 

175461553 

43.902 

16667133 

748.4 

180029010 

64.308 

8 

1 

26709 

.551 

21492 

.582 

71340013 

4.047 

461931 

073.211 

574057 

218.010 

153328107 

50101.30 

123379717 

26395.30 

19054597 

523146.2 

992809138 

5353.1400 

9 

1 

7174. 

229 

11784 

.822 

51469556 

.527 

138882 

033.857 

845470 

08.771 

606559571 

445.613 

996371464 

373.840 

36925436 

6335.649 

163670007 

3254.860 

10 

1 

14058 

.893 

13984 

.465 

19765248 

2.610 

195565 

256.251 

196606 

099.626 

276406418 

9287.230 

274943108 

3264.540 

27787751 

76073.390 

273487544 

5700.330 

11 

1 

10444 

.778 

1085. 

272 

10909339 

5.065 

117781 

4.919 

113354 

23.605 

118395987 

214.349 

123020157 

85.887 

11394563 

32395.640 

127824933 

8.993 


Table 2: Computation of Coefficient Matrix, Elements 


Computation of the model parameters, X using equation (6) 

Using coefficient matrix, A as well as the computed coefficient matrix, A elements given in table 2 and the observation 
matrix, L, the model parameters are computed as: 
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V 


"34.859077901082436" 

Cl l 


0.000225556443422 

a 2 


-0.000501234804656 



-0.000000030132587 

a 4 


0.000000003398167 

a 5 


0.000000051102750 

«6 


0.000000000004046 

a n 


-0.000000000008676 

«s 


-0.000000000000149 

\ a 9y 


. 0.000000000004302 


Substitution of the computed parameters into the interpolation model 
Having computed the model parameters, they are substituted into the geometric geoid model. Thus, 

N = 34.859077901082436 + 0.000225556443422.x - 0.000501234804656+ - 

-0.000000030132587x 2 +0.000000003398167+ 2 +0.000000051102750xy + 

+ 0.000000000004046x 2 y - 0.00000000000867xv 2 - 0.000000000000149x 3 + 

+ 0.000000000004302y 3 

Computation of the model accuracy using equation (7) 

The local geometric geoid model whose parameters have been determined was used to obtain the geoid heights of the 
points and compared with the known geoidal heights of the points to obtain the residual as shown in table 3. 


Station 

Known Geoidal 
Height, N (m) 

Model Geoidal 
Height, N (m) 

Difference in Geoid 
Height, AN 

ESO01 

34.048 

34.047 

0.001 

ESO02 

34.895 

34.948 

-0.053 

ESO03 

34.377 

34.376 

0.001 

ESO04 

34.584 

34.584 

0.000 

ESO05 

34.279 

34.281 

-0.002 

ESO06 

34.286 

34.289 

-0.003 

ESO07 

34.282 

34.282 

0.000 

ESO08 

34.375 

34.384 

-0.009 

ESO09 

34.607 

34.608 

-0.001 

ESO 10 

34.456 

34.457 

-0.001 

ESO 11 

34.175 

34.175 

0.000 


Table 3: Known and Model Geoidal Heights and Residuals 


" 0.001 
-0.053 
0.001 
0.000 


- 0.002 


V = 


- 0.003 


V T =(0.001 -0.053 0.001 0.000 -0.002 -0.003 0.000 -0.009 -0.001 -0.001 0.000 


0.000 


-0.009 


- 0.001 


- 0.001 


0.000 
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V T V = (0.001 -0.053 0.001 0.000 -0.002 -0.003 0.000 -0.009 -0.001 -0.001 0.000)x 


'o.oor 

-0.053 
0.001 
0.000 
- 0.002 
-0.003 
0.000 
-0.009 
- 0.001 
- 0.001 
v 0.000 y 


RMSE = Accuracy = ±. 


0.003 

11 


±0.0165m 


0.003 


The computed local geometric geoid model accuracy is ± 0.017m.This implies that geoid heights can be interpolated within 
the application area with accuracy of + 17mm. 


3. Conclusion 

Considering that fitting a geometric geoid surface to a set of points whose geoidal heights are known to enable the 
interpolation of the geoid heights of new points if the variables are known requires the determination of the model 
parameters using the least squares technique which is rigorous and difficult to apply. This paper has presented the step by 
step procedures to be followed when determining geometric geoid model parameters and its accuracy using the least squares 
technique. The procedures were demonstrated using a numerical example. Considering the detailed procedures and the 
numerical example, it is certain that the difficulty in the application of least squares adjustment technique in local geometric 
geoid model parameters and accuracy determination has been simplified. 
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